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Abstract. We present a simple mathematical framework for the description of the 
dynamics of glassy systems in terms of a random walk in a complex energy landscape 
pictured as a network of minima. We show how to use the tools developed for the 
study of dynamical processes on complex networks, in order to go beyond mean-field 
models that consider that all minima are connected to each other. We consider several 
possibilities for the transition rates between minima, and show that in all cases the 
existence of a glassy phase depends on a delicate interplay between the network's 
topology and the relationship between energy and degree of a minimum. Interestingly, 
the network's degree correlations and the details of the transition rates do not play 
any role in the existence (nor in the value) of the transition temperature, but have an 
impact only on more involved properties. For Glauber or Metropolis rates in particular, 
we find that the low-temperature phase can be further divided into two regions with 
different scaling properties of the average trapping time. Overall, our results rationalize 
and link the empirical findings about correlations between the energy of the minima 
and their degree, and should stimulate further investigations on this issue. 
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1. Introduction 

In the last decade, studies about the structure and dynamics of complex networks have 
blossomed, thanks in particular to the versatility of the network representation, which 
has turned out to be adequate for systems as diverse as the Internet or social networks. 
A large body of knowledge about the empirical description of networked systems has 
thus been accumulated, together with a wealth of modeling techniques; a good level of 
understanding of how dynamical processes taking place on networks depend on their 
structure has been as well reached Many network studies have been 

concerned with systems of interest in several scientific areas a priori remote from physics 
(social sciences, biology, computer science, epidemiology, ...), and they have also reached 
more traditional fields of statistical physics, such as the study of glassy systems, as we 
now describe. 

The many puzzles raised by the glass transition, and in particular the slow dynamics 
displayed by glassy systems at low temperatures, have been the subject of a large interest 
in the past decades 0, H|. One of the approaches which has led to promising insights 
consists in the description of the dynamics of a glassy system inside its configuration 
space. The energy landscape of a glassy system is typically rugged, made of many local 
minima (metastable states), whose huge number makes it difficult to reach equilibrium. 
In this framework, the energy landscape is seen as a set of basins of attractions of local 
minima ("traps"), and the system evolves through a succession of harmonic vibrations 
inside traps and jumps between minima (si, Eo|- This picture has stimulated the 
definition and study of various simplified models of dynamical evolution between traps, 
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in order to reproduce the phenomenology of glassy dynamics 
On the other hand, several studies have focused on obtaining a better understanding of 
the structure of these local minima. A way to attain this goal is to perform numerical 
simulations of small systems, at a fixed temperature, and quenching them at regular 
time intervals in order to make them reach the nearest local minimum. Information is 
then gathered on the various local minima, and on the sizes of their basins of attraction. 
Various studies have investigated, among other issues, the detailed structure of the 
potential energy landscape, the substructure of minima, and the properties of energy 



barriers between minima |18|, ll9|, [20[ . Several works have also used the information on 
the energy landscape to study a master equation for the time evolution of the probability 
to be in each minimum. The considered systems range from clusters of Lennard- Jones 

atoms to proteins or heteropolymers 0,0,0, H,£j. 

An interesting property of the modeling of the energy landscape in terms of a set 
of traps linked by energy barriers, lies in the possibility to define and study its network 
representation within the context of network theory. In this representation, each local 
minimum is associated to a node, and a link is drawn between two nodes whenever it is 
possible for the system to jump between the basins of attraction of the corresponding 
minima. The links can then be defined as weighted and directed, as jumps between 
minima are not equiprobable, and may be easier in one direction than in the other. 
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Networks of local minima of the energy landscape have thus been built and studied. 



These networks have been found to exhibit a small- world character [24[]. The number 
of links of each node (its degree) turns out to be strongly heterogeneous, possibly with 
scale-free degree distributions, which have been linked to scale-free distributions of the 
areas of the basins of attraction 25|, |26|, |27j. Complex network analysis tools have 



also been used to investigate the structure of energy landscapes of various systems 
of interest, such as Lennard- Jones atoms, proteins, or spin glasses, among others 

21 . 22, 23, 27, 28, 29, 30, 31, 32, 3^. The energy of a minimum and its degree (i.e., the 
number of other minima which can be reached from this minimum) have been shown 
to be correlated, as well as the barriers to overcome to escape from a minimum. In 
particular, a logarithmic dependence of the energy of a minimum on its degree has been 
exhibited, as well as energy barriers increasing as a (small) power of the degree of a node 

23L |25|, [271]. No systematic study of these issues has however been performed, and most 
investigations have been limited to relatively small systems because of computational 
limitations. 

Most importantly, the investigations cited above have focused on the topology of 
the network of minima, conceived as a tool to characterize the energy landscape. The 
structure of a network has however a deep impact on the properties of the dynamical 
processes which take place on it 0. It seems thus adequate to put to use the tools and 
techniques developed for the analysis of dynamical processes on networks to achieve a 
better understanding of how the energy landscape structure, represented as a network, 
affects the system performing a random walk in it, and how the onset of glassy dynamics 
can be described in this way in a general framework. In a previous paper 34J , we have 
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made a first step to fill this gap by focusing on the trap model put forward in Ref. 
In this paper, we generalize our approach to more involved transition rates between 
energy minima. We show how the heterogeneous mean-field (HMF) theory {(], [35| can 
be used in this context to highlight the connection between the topological properties 
of the network of minima and the dynamical exploration of these minima. We show in 
particular that the relationship between energy and degree of the minima is a crucial 
ingredient for the existence of a transition and the subsequent glassy phenomenology. 
Our results shed light on the empirically found relationship between the energy of a 
local minimum and its degree, and we hope that they will stimulate more systematic 
investigations on this issue. 

We have organized our paper as follows: In Section [2] we define our model of 
energy landscape dynamics as a random walk on a complex network. Different physical 
transition rates are proposed, and the corresponding numerical implementation is 
discussed. In Section [3] we present a theoretical analysis based on the heterogeneous 
mean-field approximation for dynamical processes on complex networks. This formalism 
is applied in Section HI where general analytical approximate expression are presented for 
the main quantities characterizing the glassy transition and dynamics. These expressions 
are applied to the different physical transition rates considered in Section [5j where 
checks against numerical simulations are also presented. In Section [6] we discuss the 
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Figure 1. Potential energy landscape description. Energies Ei are measured positive 
downwards. Energy gaps Sy are denned positive upwards. 

relation between energy basins and energy barriers. Finally, in Section [7] we present our 
conclusions. 

2. Random walk models on complex energy landscapes 

2.1. Definition 

We consider a network of N nodes, in which each vertex % corresponds to a minimum 
in the energy landscape, and a link is drawn between two minima % and j if the system 
can jump directly from % to j. To each node % is associated the energy —Ei of the 
corresponding minimum (energies are defined from a reference level, in such a way that 
Ei > for all i). Moreover, an energy gap Ejj is associated to the edge between 
vertices i and j, as depicted in Fig. [1] E^ is a symmetric function, such that the energy 
barrier that must be overcome to jump from vertex i to vertex j can be written as 
AEij = Ei + Sjj and, analogously, AEji = Ej + Ey. Obviously, we will have in general 
AEtj ± AEji. 

The system under investigation is pictured as a walker exploring the network 
through a biased random walk. The rate (probability per unit time) r^j to go from 
vertex i to vertex j depends a priori on the energy at vertices % and j and/or the energy 
barrier between % and j that must be overcome. The random walk model is defined in 
discrete time t as follows: 

• At time t, the walker is in vertex %. 

• It chooses at random a neighbor of i, namely j. 

• With a probability r^j, that depends on the energy Ei and/or on the energy barrier 
AE^, the walker hops to vertex j. 

• Time is updated t — > t + 1. 
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The relationship between the probabilities r^j and the energy and energy barriers 
can be of different forms. In usual unbiased random walks, r^j is a constant independent 



of both i and j [36J. As a first step to introduce a dependence on the nodes, a 
possible approach is that in which the energy barriers depend only on the local minima 
themselves, i.e. we consider Sjj = 0. For example, in the Bouchaud trap model 



considered in Ref. ll|, the probability to exit from a trap is just an Arrhenius law 



depending only on the departing trap's depth, namely 

rt^ = r e~^ 7 (1) 

where = 1/T is the inverse temperature and r is a constant that determines a global 
timescale. Other possible definitions include the Metropolis one 

r Met r0P 0Us = rQ ^ ^ e m-E^ ^ (2) 

and the Glauber rate 

Glauber _ ^0 

'^3 ~ 1 _|_ e -P{Ei-Ei) ■ 

We note that the rates considered in the Bouchaud trap model are quite different from 
the Metropolis or Glauber rates. Indeed, while the former depends only on the depth 
of the originating trap, the latter depend also on the energy of the arriving vertex. 
This translates in the fact that, in the limit of zero temperature, the dynamics of the 
Bouchaud trap model is frozen for an y E j, while Metropolis and Glauber dynamics still 



allow jumps to lower energy minima 13|. Within an even more realistic representation 
of glassy dynamics, one can also contemplate the case E^- ^ 0, allowing for the transition 
rates to depend explicitly on the energy barriers between adjacent minima. As a 
paradigm of this choice, we propose a rate of the Arrhenius form 



barriers — fiAE. 



r e~^, (4) 



which acts a straightforward generalization of the local transition rate ([Q). 

The case of rate (CD) (local trapping) was studied in a previous publication 34]. In 
the following, we will consider in turn non-local rates (j2J), (EJ) and (OS) and discuss the 
fundamental differences due to the introduction of energy barriers in the model. 

We emphasize that our model differs both from usual unbiased random walks, as 
the local energy determines the transition rates, and from mean-field trap models in 
which jumps between any pair of energy minima are a priori possible: Here, the system 
can jump only between neighboring nodes. The dynamical evolution depends therefore 
both on the network topology and on the energies associated to the nodes. 

2.2. Numerical implementation 

To implement numerically the random walk, it is convenient to resort to the techniques 
developed for general diffusion processes on complex weighted networks [37]]. The main 
advantage of this method is to avoid rejection steps, thus improving dramatically the 



computational efficiency 38|,|39|]. Therefore, at each simulation step the random walker 
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sitting at node i selects a neighbor j with probability r^j/ J2j r i->j, where the sum 
in the normalizing factor is extended to all of i's neighbors. As the walker hops on 
node j, the physical time is incremented by an interval At drawn from the exponential 
distribution P{At) = 1/At exp(- At/At), where At = h/T.3 

r is the inverse of the 

average escape rate out of node i. In this way the simulation time is disentangled from 
the physical time and the latter has no impact on the simulation efficiency No matter 
how much physical time a walker spends in a node, from the simulation time point of 
view it is always just a time step. 

The network substrates on which we will focus are scale-free networks with a 
degree distribution of the form P(k) ~ /c -7 and 2 < 7 < 3. We will generate them 



with the Uncorrelated Configuration Model (UCM) [40J, that allows us to tune the 
degree distribution to the desired form and prevents the formation of degree-degree 
correlations. Networks are therefore generated as follows: A number of stubs (or semi- 
links) extracted from the desired final degree distribution is assigned to each node. 
Stubs are then randomly paired to form links between nodes, with the prescription that 
multiple links as well as self-loops must be avoided. A minimal degree m is fixed a 



priori. To avoid spurious effects due to the possible presence of tree-like structures 41 



it is convenient to adopt m > 2. We will choose m = 4 in all of our simulations. So 



far the algorithm coincides with the one of the Configuration Model 42|, but the UCM 



introduces moreover a cutoff to the degree distribution, k c = N 1 ^ 2 , which avoids the 
formation of degree correlations by limiting the size of the hubs 40[ . 

3. Heterogeneous mean-field theory 

In order to gain analytical understanding of the role of the different transition rates 
in the corresponding glas sy dynamics, we apply a standard heterogeneous mean-field 



(HMF) formalism p, |35[ . The basic tenet of HMF is the assumption that all the 
dynamical properties of a vertex depend only on its degree. Vertices are thus grouped 
into classes according to their degree, and vertices with the same degree are treated 
as equivalent. This approximation is consistent with previous findings that have 
uncovered the correlations between the energy of a local minimum and the degree of the 



corresponding node in the network [25[. We therefore make the assumption that there 
exists a relationship Ei = h(ki) where the function h(k) is a characteristic of the model. 
This also means that the distributions of energies p{E) of the system's landscape, and 
the degree distribution P(k) of the corresponding network are linked through h. In the 
same spirit, we make the further assumption that the energy gap between minima i and 
j depends only on the degrees of i and j, i.e., that it can be written as = o~(ki, kj), 
where a(k, k') is a symmetric function of k and k' . 

Under the HMF approximation the dynamics will thus focus on the transitions 
between different degree classes. The rate to go from a vertex A; to a vertex k! can be 
written as 

W kk , = P(k'\k)r(k -)• k'). (5) 



Complex networks and glassy dynamics: walks in the energy landscape 



7 



The function P(k'\k), defined as the conditional probability that a vertex of degree k 



is connected with another vertex of degree k' |43j, takes into account the topological 
features of the network, by gauging the probability of selecting a vertex k! as neighbor 
of k. The function r(k — > k') measures the rate of jumping from a vertex of degree k to 
a vertex of degree k' (given that they are connected by an edge) , and depends on k and 
k! through the rates r^j and the functions h and a. Obviously, the rate r(k — > k') is 
not in general a symmetric function of k and k'. It is worth noting that, apart from a 
normalization, Equation ((5]) is simply the so-called weighted propagator describing the 



probability that a node in class k interacts with a node in class k' |37j]. We also note 
that the rates rik — > k') depend on the inverse temperature /3 through the microscopic 
rates r^ ^j. 

4. General HMF formalism 

In this section, we apply the HMF theory to compute different quantities relevant for 
the characterization of the dynamics of a random walk in a complex energy landscape 
represented in terms of a network of minima. 

4-1. Occupation probability 

The description of a random walk dynamics starts from the occupation probability 
P(k,t w ), defined as the probability for the walker to be in any node of degree k at a 
time t w . Its time evolution can be easily represented in terms of a master equation of 
the form 

dP{ ^ tw) = P(k, = W w P(k, t w ) + £ W k 'kP(k', t w ) . (6) 

™ k' k' 

Upon describing the state at time t w with the row vector P(t w ) = 
{P(l, t w ), P(2, £„,),.., P(k c , t w )}, where k c is the cutoff or largest degree in the network, 
Equation ([6]) can be rewritten in vector form as 

P(t w ) = -P(t w )L , (7) 

where the matrix L, with elements 

L k 'k = (5k'k J2 W ki - W k'k^j , (8) 

is a generalization of a Laplacian matrix to the case of directed weighted graphs. The 
matrix elements satisfy 

Lk'k' = X! Lk'k , (9) 
which ensures conservation of probability and states that the columns of L are not 



linearly independent. The real part of every eigenvalue of L is non- negative 44(. As a 
consequence, all solutions of Eq. ([7]), which can be formally written as 

P(0 = P(to)e" L(t -' o) , (10) 
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are stable according to Lyapunov criteria. In particular, since det L = 0, L always has 
the eigenvalue 0, which corresponds to a constant solution of the problem. At this point 



one can proceed in close analogy with discrete-time regular Markov chains [45|. By 
making the assumption that the matrix W k > k is non-negative and irreducible (indeed it 
is for every choice of rikl — > k) in the following), we can prove that the eigenvalue of 
L has algebraic multiplicity 1. Hence, the stationary solution of Eq. ([7]) is unique. 

4-1.1. Steady state In order to calculate the steady solution P°° in the limit t w — > 00, 
one can impose P(t w ) = 0. This leads to the condition 

P°°L = 0, (11) 

so that we are left with the task of finding the left nullspace of L. Eq. (1TT1) is a 
homogeneous system of algebraic linear equations. It admits non-trivial solutions since 
det(L) = 0. In our case, the solution to (TTTT) can be easily found by imposing the 
detailed balance condition. Namely, writing Eq. (TTTT) as 

]T [-W w P°°{k) + W k , k P°°(k')} = 0, (12) 

k' 

we can obtain a solution by imposing that the terms inside the summation in Eq. (1121) 
cancel individually, that is 

W kk/ P°°(k) = Wk /k P°°(k'), Wk, k'. (13) 

Substituting the form of W kk >, we obtain 

P°°(k) _ W k >k _ P(k\k')r(k' -» k) _ kP(k) r(k' -> k) 
P°°(k') ~ Wk~k~> ~ P{k'\k)r(k -> V) ~ k'P(k') r{k V) ' ' 
where in the last step we have used the degree detailed balance condition kP{k)P{k'\k) = 
k'P(k')P(k\k') which simply expresses that the number of edges from a node of degree 
A; to a node of degree k! is equal to the number of edges from a node of degree k' to a 



node of degree k [46|. From Eq. (1141) . we see that its right-hand-side must be expressible 
as a simple ratio of a function of k over a function of kl . A general way to obtain this 
is to impose a coarse-grained rate r(k' — > k) taking the general form 

r(k' -+k) = f(k')g(k)s(k',k). (15) 

In other words, we assume that the rate r{k' — > k) can be written as the product of a 
function of k', a function of k, and a symmetric function s(k', k) = s(k, k') (where k and 
k' need not be separable). We will see later that all the rates r^j defined in Sec. 12.11 
(traps, Glauber, Metropolis, and energy barriers) can be written in such a form. The 
stationary solution is then given by 

P°°(k) = l k P(k)g(k)/f(k) (16) 

where Z is a normalizing constant determined by the condition J2 k P°°{k) = 1- Such 
a solution is unique, as proven above. Interestingly, the symmetric function s(k',k) 
does not enter the steady solution, although it will play a role in affecting the transient 
behavior, as we will see in the next sections. 
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4-1.2. Glassy phase The steady state solution found above for the occupation 
probability is defined if and only if the normalization constant 

Z = J2kP(k)g(k)/f(k). (17) 

k 

is finite. When this condition is met, the random walker reaches an equilibrium state 
with a distribution P eq = P°°. On the other hand, whenever such condition is not met, 
the random walker is unable to reach a steady state, i.e. the steady solution to the 
rate equation does not correspond to any physical steady state in equilibrium P eq . We 



identify this region of the phase space with the glass phase for our random walker |H 

The functions / and g depend on the temperature, on the precise dynamics chosen 
(traps, Glauber, Metropolis), and encode the relationship h between energy and degree 
of the minima. The degree distribution moreover enters explicitly the expression Z. 
As the various parameters of the model are changed, it is thus a priori possible to go 
from one phase in which Z is finite to one in which Z diverges. In a physical system 
in particular, the control parameter is usually the temperature, while the topology of 
the network of minima and the function h are given. It is then clear from Eq. ( jTTf) 
that the presence or absence of a finite glass transition temperature f3 c , such that Z 
becomes infinite for (3 > f3 c , depends on the interplay between the topology of the 
landscape network (as determined by P(k)) and the functions / and g. Interestingly, 
at this mean-field level, the existence of a transition does not depend on the network 
degree correlations, since the conditional probability P(k'\k) do not enter Eq. ( ITTj) . 

Let us consider for instance a network of minima with a heavy-tailed degree 
distribution such as P{k) ~ A; -7 . A transition between a finite and an infinite Z 
can be observed if and only if g(k)/f(k) behaves at large k in the form ~ k a where the 
exponent a depends on the temperature, and can take values smaller or larger than 7 — 2 
depending on the temperature. Another example is given by a stretched exponential 
form for P{k), P{k) ~ e~ 6fca , in which case a transition is observed if and only if 
g(k)/ f(k) is of the form e b ' fca , with b' a function of the temperature (the transition is 
then given by b'((3 c ) = b). 

4-1.3. Glassy dynamics In any finite system, unless the product function g(k)/f(k) 
exhibits some sort of singularity, the normalization constant Eq. ( jTTl) is finite and the 
steady state distribution P°°(k) exists, the occupation probability P{k,t w ) converging 
to it after an equilibration time, i.e. 

lim P(k,t w ) = P°°(k). (18) 

The corresponding thermalization of the occupation probability occurs in a way 
depending on the function h. Shallow energy minima are indeed explored first, while 
deep traps (large E) are visited at larger times [ill, 0] • If /i is a growing function of 



k, as indeed found empirically [25|, small degree nodes correspond to shallow minima, 
and deeper minima are associated to larger nodes. The evolution of P(k,t w ) takes 
then place in a hierarchical fashion: The small degree region equilibrates first, and 
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progressive equilibration of larger degree regions takes place at larger times. In this 
respect, we obtain a strong difference between the biased random walk that the glassy 
system experiences and usual diffusion processes corresponding to unbiased random 
walks, which visit first large degree vertices and then cascade down towards small degree 
nodes 36, 47, fjj], in the present case we observe an "inverse cascade process" from small 



vertices to hubs. 

We have found in |34j that, in the case of a random walk among traps, this 
hierarchical thermalization is summarized in a scaling form for P(k,t w ), which can 
be written as 

p(k, t w ) = k w (t w y l F ij^pjj > (is) 

where k w {t w ) represents the maximum degree of the vertices equilibrated up to time t w , 
and F{x) interpolates between P°°(x) at small x and the short time form of P(k,t w ) 
which is proportional to kP(k). We will see in the next section that a similar scaling 
is obeyed for other transition rates. In general, for the glassy dynamics, the functional 
form of k w (t w ) can moreover be obtained through the following argument: The total 
time t w can be written as the sum of the trapping times spent in the vertices that have 
been visited since the beginning of the dynamics. Trapping times increase with the 
depth of the minimum, hence with the degree (we are still considering the case of an 
increasing function h(k)), and, in the glassy phase, the consequence is that the sum of 
trapping times is dominated by the vertex with the largest degree visited up to that 
point, namely k w . Moreover, the average trapping time r k at a given vertex k can be 
estimated as the inverse of the average escape rate r k from that vertex: 

- = r k = ]T W kk , = ]T P(k'\k)r(k = £ P(k'\ k)f(k)g(k')s(k, k').(20) 

Tk k' k' k' 

We can therefore estimate k w (t w ), the typical degree up to which nodes are 
"equilibrated" at time t w , by approximating r kw ~ t w , and solving the equation 

11 

tw ~ f(k w )Ek'P(k'\k w )g(k')s(k',k w ) (21) 
to obtain k w as a function of t w . Note that the result depends here on the function 
s(k, k') and not only on /, g, h and P(k). 

4-2. Average escape time 

The properties of the system can be further quantified by measuring the average time 
tpsr(t m ) required by the random walker to escape from the vertex it occupies at time t w 



34J. For small waiting times t w , t esc increases as a result of the transient equilibration 
of P(k,t). For large t w , such that P(k,t w ) is close enough to the equilibrium P°°, t esc 
can be calculated instead as the average 

t esc (t w oo) = J2P°°(k)r k = ^Y, kP ( k Mk)/f(k)]r k (22) 
k ^ k 
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where T\~ = 1/r^ is the inverse of the equilibrium escape rate, cf Eq. (|2"U|) . yielding 
i d , ^ 1 V kP(k)g(k) 

eSc[w fW 2 Eit' P(A:1%( W, fc) • ^ 

Most interestingly, the explicit form of the average escape time t esc depends explicitely 
on the symmetric function s(k, k') as well as on the network degree correlations, as 
expressed by the conditional probability P{k'\k). 

4-3. Average rest time 

Let us go back to the issue of the existence of a glass transition in the model. We 
first recall the phenomenology of the fully connected trap model, with transition rates 
rj_>j = roe~ l3Ei /N for any i and?, where the energies Ei are random numbers extracted 
from a distribution p(E) [HI, Q. As all traps are connected with each other, all traps 
are equiprobable after a jump, so that the probability for the system to be in a trap of 
depth E is simply p(E), and the average rest time spent in a trap is (r) = / p(E)e l3E dE. 
A transition between a high temperature phase and a glassy one is thus obtained if and 
only if, when f3 increases, (r) is finite at small /3 and diverges at a finite j3 c . Such a 
phenomenology is obtained if and only if p(E) is of the form exp(— /3 C E) at large E (else 
the transition temperature is either or oo), and the transition temperature is then 
T c = 1/& |ll 



In the present case of a nework of minima, the average rest time that the walker 
spends in a minimum is 

= (24) 

where the symbol (...)h refers to the average performed over the measure Ph{k), which 
represents the probability that the walker is in any vertex of degree k after a hop. Note 
that we disregard here the physical time, which is the sum of times spent in the various 
minima, but consider only the number of hops between minima. In the case of the traps 
model, Ph is simply given by the probability to be in a node of degree k after a hop in 



a random walk, i.e. by kP(k)/(k) [34|, since the transition rates do not depend on the 
arrival node. In a non-local trapping model instead, we need to write a master equation 
of the form 

Ph(k) = -P h (k) + J2W k 'kPh(k'), (25) 
k' 

where the matrix Wk'k = Wk'k/ J2i Wyi = Wk'k/fk' is now stochastic and the derivative 
is intended with respect to the number of hops. In the long time limit, we impose 
Ph(k) = and calculate Ph{k) as we did for P°°(k), imposing the detailed balance 
condition, and obtaining 

P h (k) = ^kP(k)[g(k)/f(k)]r k , (26) 
where X is a normalization factor, given by 

1 = E E kP(k)P{l\k)g(k)g(t)s{k, I). (27) 
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We finally obtain for the average (r): 

0-> = £^(*)A-* = f' (28) 

k x 

where Z is the normalization factor of P°°(k) defined in Eq. ( fl7l) . As for the average 
escape time, the average rest time (r) thus depends on all the parameters of the system, 
including the network's degree correlations and the symmetric function s. 



5. Application to physical transition rates 



In this section, we apply the general HMF results obtained in the previous section to 
physical transition probabilities between local minima given by the trap model, Glauber, 
Metropolis and barrier-mediated rates. We will focus for definiteness on scale-free 
networks characterized by a power-law degree distribution P(k) ~ fc~ 7 with 2 < 7 < 3, 



23 



25 



which turns out to be the interval reported in the literature 

Let us first consider the explicit form of the transition rates in each case, to show 
that they can be cast in the form given by Eq. ({15]) . In the case of the trap model, the 
rate to jump from a vertex k to a vertex k' is simply r{k — > k') 



r e 



r e 



-Ph(k) 



where we recall that h(k) gives the depth of a node of degree k: it depends only on 
the degree of the starting node, and not on the node reached after the jump. We can 
therefore use 



(Trap model) f(k) = e 



-Ph(k) 



g(k) = 1, s(k,k') = r . 



The Glauber rate can be written as 



r(k — > k') 



oPh(k') 



r 



e l3h{k) _|_ e Ph(k') 



leading to 



(Glauber) f(k) = 1, g(k) = e ph{k \ 



s(k,k') 



r - 



■PHk) _|_ e /3h(k') ■ 



The Metropolis transition, on its turn, reads 

r(k k') = r min [l, e PW)-h{k)) _ 

Since, for positive a, min(l,6/a) = min(a, b)/a, we can choose 
(Metropolis) f(k) = e~ pm , g{k) = 1, s(k,k') 



(29) 



(30) 



(31) 



(32) 



ro min( 



Finally, in the presence of energy barriers, the transition rate reads 



r(k — > k') 



r e 



-/3{h(k)+a(k,k')) 



where a(k, k') is a symmetric function of its arguments, so that we can use 



(Barriers) f(k) = e 



-/3h(k) 



g(k) = l, s(k,k , )=r e-^ k '\ 



(34) 



(35) 
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5.1. Steady state and the glass transition temperature 

Interestingly, for all the transition rates considered above, the product of the functions 
1/ / and g, which controls the existence of the steady state solution of the occupation 
probability, takes the form 

g{k)/f{k) = e^ k \ (36) 

The normalization constant Z defined in Eq. ( 1T7|) can thus be written as 

Z = J2kP(ky h{k \ (37) 

k 

For a power-law degree distribution P(k) ~ /c -7 , a finite glass transition temperature is 
then obtained if and only if h is of the form 

h(k) = E \og(k) , (38) 

which is precisely what has been found, in conjunction with a scale-free degree 



distribution, in Ref. [25]. Z is then indeed given by a sum of terms of the form 



k 1 ~< + P E ° , which converges if and only if 

(3E - 7 < -2 . (39) 
In other words, a transition between a high temperature phase in which P°°(k) exists 



and a low temperature glassy phase is obtained at the critical temperature 34] 



Quite noticeably, the existence of a transition at a finite temperature, as well as 
the value of this temperature, does not depend on the form of the transition rates 
between the local minima, but only on the existence of a particular interplay between 
the topology of the network of minima and the relationship between energy and degree 
in this network, as determined by the function h. We emphasize that this result is also 
independent of the network degree correlations P(k'\k), as already noted in the previous 
section. 

5.2. The steady state and finite size effects 

Let us focus on the case of a scale-free network of minima, with P(k) oc /c~ 7 and 
Ei = Eolog(ki). For any of the rates discussed above, the steady state measure, when 
it exists, is given by 

U-7+/3E0 

^m = cM77^y for7 -^ >2 ' (41) 

where ( is the Riemann ( function. A plot of P°°(k) as a function of k and 7 — f3E 
is given in Fig. |2l while data from simulations are reported in Fig. [3] for the evolution 
of P(k,t w ) under Glauber dynamics. Above the transition (7 — (3E > 2), low-fc states 
(i.e., shallow minima) are more probable. As the temperature decreases, P°°(k) becomes 
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Figure 2. Equilibrium probability distribution P°°(k) for the random walker to be in 
any node of degree k. For 7 — /3Eq = 2 the system undergoes a transition to a glassy 
state. 



less and less peaked at low values of k, and large-fc states, which correspond to lower 
energies, become more and more probable. 

In any finite system, the sum defining Z, is finite at any temperature as the degree 
distribution has a cut-off at a finite k c : 

Z = Y.kP{k)g{k)/f{k). (42) 

k 

For instance, for P(k) oc fc -7 , and with h(k) = Eo\og(k), the sum 

£ tf-T+PEo = H (-1+1-PE ) (43) 
k=l 

is analytic in 7 — /3E = 2 for any finite k c . Here is the Harmonic Number of order 
a, which tends to C(°0 f° r k c — > 00. 

The probability P°°(k) is thus well defined for every 7 — f3E and for any finite 
system. In particular, performing a continuous degree approximation in Eq. ( 143|) . we 
can obtain an estimate of the network size dependence of Z as 

Z ~ [ kc k l ^ +f3Eo dk ~ const + k 2 ~^ Eo . (44) 



For < (7 — 2)/E , Z tends to a constant as the network size (and thus k c ) increases. 
On the other hand, for > (7 - 2)/E , Z diverges as k^ +l3Eo , i.e., as N^' +f}Eo ^ 2 in 
uncorrelated scale-free networks, which obey k c ~ N l l 2 . 

5.3. Glassy dynamics 

At low temperatures, even for a finite system, the evolution of P(k; t w ) towards P°°(k) 
is slow, as displayed in Fig. [31 and an ageing regime takes place, in which the function 
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Figure 3. Evolution towards equilibrium of the probability distribution P(k;t w ) 
for Glauber dynamics. The distribution measured after a small waiting time t w is 
determined by a usual unbiased random walk behavior, i.e P(k;t w ) ~ kP(k), while at 
larger times it relaxes to the equilibrium P°°(k). The relaxation towards equilibrium 
starts from small degree nodes. Data refer to UCM networks with N = 10 6 and 7 = 3.0 
for (3 = 2 (Eq = 1): with these parameter values, for small times P(k; t w ) ~ fc~ 2 , while 
at large enough times P(k; t w ) ~ P°°(k) = const. 



P(k,t w ) obeys the scaling form 

P^k^tyj} k w {t w ^j J~ 



k 



(45) 



where the characteristic dee; be estimated from Eq. (121]) . In order to 

simplify its computation, we will consider an uncorrelated network of minima, such 
that P(k'\k) = k'P(k')/{k), and we will work in the continuous degree approximation, 
using the normalized form P{k) = (7 — l)m 7_1 /c~ 7 , where m is the minimum degree 
present in the network. 

In the case of the Glauber dynamics, the escape rate can be expressed, within the 
above approximations, as 



1 



k'P(k')r 



7-2 7-2 
1+ ' 



D /3h(k') 







e Ph{k) _^ e Ph(k') 

(-)' 



dk' = 



7-1 



(I 



(k) 



z 



dz 



(46) 



48 



where we have used the relation h{k) = E \nk and where T[a, b, c, z] is the Gauss 
Hypergeometric function. Using the asymptotic expansion for z — > 
that the leading behavior for large k yields 

' k~^ Eo 0> C 
k-^ 0<0 C ' 



we obtain 



(47) 
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,3 



Figure 4. Data collapse for the time evolution of the occupation probability P(k; t w ) 
at different temperatures (Glauber dynamics). Data refer to UCM networks with 
N = 10 6 and 7 = 2.5, so that (3 C = 7 — 2 = 0.5 (where we have taken Eq = 1). The top 
panel presents data for /3 < (3 C while both middle and bottom panels concern the low 
temperature case /3 > /3 C . Accordingly, for the top panel we use kyj ~ tw^ ^ for 
the rescaling, while for the central and the bottom ones it holds k w ~ <i = (see 
Eq. (|49p ). The curves corresponding to different t w collapse well under this rescaling. 
Note that we use rather small values of t w , because the equilibration time, defined by 
k w ~ k c ~ N 1 / 2 , is t eq ~ ~ 32 for /3 > (3 C and t eq ~ ~ 6 for (3 < f3 c . Each 
curve is obtained by averaging over 3 x 10 6 simulation runs. 



which leads to 

1 / p > ft 

T * = ^~{^» ^<ft. • (48) 

From here, using the relation Tk w ~ t w , we obtain 

In Fig. H] we check the validity of Eqs. (145 p and (I49p by performing a data collapse 
analysis for different values of t w . The curves obtained for different t w collapse indeed 
as predicted. 

In the case of the Metropolis transition rates, a similar analysis yields 

Metropolis 1 [ U~PE P h -/3 c E \ (z^ 

K Wc - mo { Wc J' (50) 

leading to the same asymptotic behavior as in Eq. ( j47|) . and therefore to the same scaling 
picture as for the Glauber rate. 

As pointed out for the case of the traps model 34|, however, in finite systems the 
scaling relations in Eq. ( H9|) hold only as far as k w {t w ) < k c , i.e. it exists an equilibration 
time t eq , obtained by inverting Eq. P9|) . above which the system has completely relaxed 
and Eq. (T4"5"l) is no more valid. Finally, it is worth stressing that, while for large 
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temperatures the scaling exponent relating k w and t w depends on the temperature, in the 
low temperature phase it becomes independent, being just proportional to the transition 
temperature. We note that this saturation of the exponent at 1/(/3 c Eq) is very different 
from the phenomenology obtained in the trap model 34j, for which k w ~ t\l^ E °\ 
An immediate consequence is that the equilibration time strongly depends on /3, as 
teq ~ k^ E °\ for a system described by traps, but is given by t eq ~ k^ cE °' <C k^ E °^ for 
any /3 > j3 c for Glauber and Metropolis rates. 

Contrarily to results for the glass transition temperature T c and the steady state, 
the glassy dynamics for barrier-mediated rates does not yield the same results as for 
Glauber and Metropolis rates, since does depend on the symmetric function a(k, k'). 
In particular, we need here to choose a functional form for a. We propose to use 

a{k,k') = a (k^ + k^), (51) 

which will be justified in Section EJ In this case we obtain 

farriers ^ ^-/HSog-^k^ ( 52 ) 

In order to keep an interesting phenomenology, the constant o"o cannot be chosen 
arbitrarily. If o"o were independent of the system size, the escape rate would be 
dominated by the exponential at all temperatures. This behavior would reflect the 
fact that in this case the rate is suppressed in transitions involving nodes of large k, 
eventually generating unphysically large barriers at large k c . To prevent the system 
from building up infinite barriers, one can impose 

E lnk c ota k%, (53) 

such that the maximum barrier is always comparable to the lowest energy minimum 
and neither term dominates the other. As a consequence, we take 

a = € E — , (54) 

where e is now constant and size independent. Contrarily to the previous cases, k w (t w j 
is hard to determine, as no explicit inversion of Eqs. (|21J) and fl52|) can be provided for 
the range of parameters of interest in our study. A numerical evaluation of k w (t w ), is 
reported in Figure [5j The maximum degree of equilibrated vertices k w has an initial 
power law increase in time, which is reminiscent of the local trapping model. However, 
as larger degree nodes are equilibrated, exponential barriers come into play and the 
hierarchical thermalization becomes logarithmic in time. 



5.4- Average escape time 

The average escape time, defined as the average time required by the system to escape 
from the vertex it occupies, can be computed in the long time limit from Eq. (122j) . as a 
function of the average trapping time in vertices of degree k. From the asymptotic 
expansions of in Eq. (jUJ), valid for Glauber and Metropolis dynamics, evaluation of 
Eq. ( 1221) allows us to observe that, whenever a finite T c exists, t esc (t w — > oo) diverges 




Figure 5. Maximum degree of equilibrated nodes up to time t w for barrier- mediated 
dynamics. Curves are obtained from the numerical inversion of Eq. (|21[) . Values of 
parameters are: 7 = 2.5, (3 C — 0.5, Eq = e = 1, k c = 10 3 , \i = 0.5. 
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Figure 6. Rescaled average escape times (Metropolis dynamics). Data for UCM 
network with 7 = 3.0, so that j3 c = 1.0 (Eq = 1). The scaling forms of Eq. ([55)1 produce 
a collapse of the curves concerning different systems sizes in the three regimes of high, 
intermediate and low temperature (top, center and bottom panels, respectively). The 
slightly worse collapse obtained for f3 = 0.75 may be due to logarithmic corrections as 
/? is close to both j3 c and 2/3 c . Each point is averaged over 400 simulation runs (20 
runs on each of 20 network realizations). 



at 2T C in an infinite system, as was already observed in the case of local trapping 34 
It is noticeable that the same divergence temperature is obtained, as Eq. (I22j) a priori 
involves the network's degree correlations and the function s. Within the continuous 
degree approximation, the divergence of the escape time with the system size follows 
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Figure 7. Average rest time for the Glauber dynamics in a scale-free uncorrelated 
network with 7 = 2.75, as a function of the inverse temperature, for different system 
sizes. Data are obtained by numerical computation of Eq. (|28|) (with Eq = 1). Note 
the exponential increase with /3 for (3 C < (3 < 2f3 c , which saturates for f3 > 2f3 c as 
predicted by Eq. (f6Tl) . 



the scaling laws 

f k&>* (3 > & 

tH = tesc(t w -+ OO) ~ J P°°(k)T k ~} kW-MEo P c / 2 <P<P C .(55) 

[ const. p < p c / 2 

As noted in the previous paragraph for the equilibration time, we note that the scaling 
for /? > f3 c differs from the form k^ E ° encountered for local trapping 34j. Figure[6]reports 
simulation data that confirm the validity of Eq. ( 1551) . As the temperature is lowered, the 
initial transient becomes longer, but for large enough times t w the asymptotic behavior 
predicted in Eq. ( I55j) is reached, as made clear from the collapse of curves concerning 
different system sizes. 

In the case of barrier-mediated dynamics, Tk depends on the symmetric function 
a(k, k'), as expressed in Eq. ( |52|) . namely 

r k ~ k^ Eo e^ ok ". (56) 
Proceeding as above we obtain 

tesc-l fcf+^^o /3 c /( 2 + e)</3</3 c , (57) 
[ const. (3 < p c / (2 + e) 

where we recall that e does not depend on the system size. 

5.5. Average rest time 



The HMF expression for the asymptotic average rest time, defined as the average time 
spent by the system in a minimum, is given by Eq. (128j) . namely (r) = Z/T, where the 
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Figure 8. Average rest time for the Glauber dynamics in a scale-free uncorrelatcd 
network with 7 = 2.5, as a function of the degree distribution cut-off and for different 
temperatures. Data are obtained by numerical computation of Eq. (|28[) . (r) grows as 
a power-law of fc c , with an exponent that grows as j3 increases (going from bottom to 
top in the figure). The thick-gray line corresponds to (3 = /3 C , while the thick-black line 
corresponds to /? = 2f3 c . For larger values of /3, the power-law behavior corresponds 
to the predicted (r^) ~ k^ cEa ~ kj~ 2 , which no more depends on /3. The kj~ 2 curve 
is reported as a dashed line for reference. Values of /? represented here are comprised 
between 0.25 and 3. 

quantities Z and X, for uncorrelated scale-free networks and a degree-energy relation 
h(k) = E \n(k), take the form, in the continuous degree approximation, 

Z ~ k l ~ 1+l3Eo dk ~ const + kf-^ )E \ (58) 

X ~ J" dk J kc dk' k^gifyk'^gik'^ik^'). (59) 

Let us first recall the case of the local trap model. Both g(k) and s(k, k') are then 
constants, so that X ~ (A;) 2 = const. Thus, the average rest time behaves as Z: it is 
finite for (3 < f3 c , and diverges with the system size as k^P~ /3c ^ E ° for (3 > (3 C , signaling the 
emergence of the glassy regime at low temperatures. 

In the case of Glauber and Metropolis dynamics (both leading to the same results), 
the situation is more involved, since the product g(k)g(k')s(k, k') entering X is not 
constant. In fact, X diverges with k c for /3Eq > 2(7 — 2), that is, at a lower temperature 
given by f3' c = 2f3 c . The interplay of these two temperatures determines the behavior 
of the system for finite sizes within the low temperature phase. In particular, for the 
Glauber dynamics with g(k) = and s(k, k') = ro/[e^ h ^ + e^ fe )], we have 

X~ const + k { /- 2 ^ Eo . (60) 

Upon considering lower values of (3, (r) first encounters the divergence of Z at (3 C , which 
is then partially regularized by the divergence of X at 2/3 c . From these results, we obtain 
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the emergence of three scaling regimes for the behavior of (r) as a function of the system 
size: 

[ k^ E ° (3>2(3 C 
(t)°° = (r)(t w oo) ~ kf-^ E ° (3 C <(3<2(3 C . (61) 

I const. f3 < f3 c 




Figure 9. Average rest time as a function of j3 for different transition rates, for a 
random walker on UCM networks with 7 = 2.75 (Eo = 1) and TV = 10 6 . Discrete 
points: simulation results; continuous lines: theoretical predictions from (r) = Z/X, 
based on simulation parameters. The Glauber and Metropolis transition rates induce 
two changes of behavior at j3 c and at 2/3 c , the first being a steep increase of the average 
rest time (r) and the second a smoothing/saturation of this increase. No saturation 
of (r) is instead observed when barriers are present. The agreement with theoretical 
predictions is remarkable, thus corroborating the validity of the HMF assumptions. 
Moderate deviations are found only in the case of barriers, where exponential growth 
is expected to add greater fluctuations. In the inset, the difference between the two 
behaviors is more evident thanks to a different scale of the plot. Note that, since in 
the simulations Eq = 1, the high temperature limit of the rest time is different for the 
Glauber and Metropolis dynamics, being t(/3 = 0) = 2 and t(/3 = 0) = 1 respectively. 
For the case of barriers we have chosen Co = 10 _1 . Each point is obtained by averaging 
the rest times corresponding to the first 10 6 hops of the random walker in each of 10 
network realizations. 

The direct numerical computation of Eq. f )28|) is shown in Figs. [7] and El showing 
the validity of this analysis. In particular, the exponential increase of (r) with j3 in 
the intermediate temperature range (3 C < (3 < 2(3 C is clearly apparent in Fig. [7J and 
Fig. [8] confirms that, for (3 > 2j3 c , the exponent in the scaling law for the system size 
k c does not depend on the temperature. While the temperature /3 C signals the onset of 
the low temperature phase with glassy dynamics for all considered transition rates, for 
Glauber/Metropolis dynamics the low temperature phase can be further divided into 
two regions that correspond to different behaviors of the timescales with the system 
size. 
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Figure 10. Average rest time for Metropolis dynamics on uncorrelated scale-free 
networks with 7 = 3.0 (Eq = 1, [3 C = 1.0). Data for different system sizes 
collapse well when rescaled according to the theoretical values of Eq. (foTj). While 
the agreement is excellent both for high and low temperature (top and bottom 
panels, respectively), logarithmic corrections are probably present for the regime of 
intermediate temperatures fi c < f3 < 2f3 c (central panel). In each simulation run the 
rest interval starting before t w and ending after t w is considered, and each point in the 
Figure is averaged over 400 simulation runs (20 runs on each of 20 network instances). 



Figures |9] and [10] moreover show the result of numerical simulations of random 
walkers on scale-free networks for Glauber and Metropolis dynamics as well as in the 
case of barriers, globally confirming the above discussed picture. 

Dynamics in the presence of barriers do not yield the same phenomenology as 
Glauber and Metropolis rate. In this case, we have g(k) = 1 and s(k, k') = r e~ /3cT( - k ' k '\ 



Selecting a(k, k') = aoik^ + /c /At ), as in Section [5731 we are led to 

1 2 



1 



k 

° dkk^e-^ ^ 



(62) 



As for the escape time t esc , upon choosing a size independent, the rest time (r) will 
be diverging exponentially with k c at every temperature. By introducing the size 
dependence as in Eq. (1541) . instead, one can see that the X integral converges to a 
constant for large k c so that one is left with 

f k (f3-/3 c )E g o 

(t) ~ c o o ■ (63) 

x ' \ const. p < fi c y ' 

We therefore obtain the same picture as in the case of traps, with an exponential increase 
of (r) as (3 increases, as confirmed by numerical simulations in Fig. [9j 
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6. Energy basins and energy barriers 



Inspired by analogies with systems governed by the Arrhenius law, we have introduced 
a transition rate that takes into account the energy barriers between states. Within 
the heterogeneous mean-field approximation, in which all variables depend only on the 
degree of the vertices, and choosing E{ = Eq ln(fci), the transition rate we have considered 
becomes 

r(k k') = k-^ E ° e -^ k ' k '\ (64) 

where a(k, k') is a symmetric function of the degrees of the two nodes. This model 
represents in essence an extension of the local trap model, where non-locality enters 
only in the form of symmetric energy gaps E^- (see Fig. [[]). The steady state has exactly 
the same form as the ones discussed so far, which incidentally is the same as for the local 
trap model. As shown in the previous Sections, the presence of barriers affect transient 
relaxation phenomena, but not the steady state. 

A different question is whether one can be more specific about the realistic 
functional form of the coarse-grained function a(k, k'). In the previous section we have 
already introduced a definition of a(k,k'). Here we provide the rationale behind that 
choice. 

Numerical simulations of the energy-landscape network of Lennard- Jones clusters 
show that the average barrier to escape from state k follow the power law AE^ ~ k^, 



with fi > |23j. In our model, such average can be computed as 

AE k = Y^P(h\k) [E \nk + a(h,k)} . (65) 

h 

For simplicity we focus on uncorrelated networks, as simulations indeed show weak 
degree correlations. Under this assumption, the first term of the sum on the right-hand 
side of Eq. fl6B]l will contribute as a logarithm of k and the power-law behavior of AEk 
is possible whenever a(h, k) ~ /c M , which leads to consider the form proposed in previous 
sections, 

a(k, k') = do (F + k'^) , (66) 

where <To has the dimensions of an energy (a discussion about the possible values of 
cr is given in Section [573]) . More complicated functional forms can also be proposed, 
for example accounting for barriers of different signs, as long as they retain the same 
power-law behavior of Eq. ( 1661) in the large k limit. It is interesting to notice that 
AEk ~ k^ 1 implies that the average escape rate e~ l3AEk has the form of a stretched 
exponential ~ exp(— Pk^), if we neglect the logarithmic correction. 



7. Conclusions 



In this paper, we have presented a simple mathematical framework for the description of 
the dynamics of glassy systems in terms of a random walk in a complex energy landscape. 
We have shown how to incorporate into this picture the network representation of this 
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landscape, put forward and studied by several authors [25|, |26|, |27J, |29|, |30|, |3l|, |32|, |33 



in order to go beyond simple mean-field models of random walks between traps that 
are all connected to each other. While our previous work had focused on the case of 
a landscape consisting of traps connected by a network 34j, we have here generalized 
our study to more involved and realistic transition rates between minima, including 
Glauber or Metropolis rates, and the possibility of energy barriers between minima. 
We have shown how the interplay between the topology of the network of minima 
and the relationship between the energy and the degree of a minimum may determine 
a rich phenomenology, with the existence of two phases and of glassy dynamics at 
low temperature. Interestingly, the existence of these phases, and the transition 
temperature, do not depend on the network's degree correlations nor on the precise form 
of the transition rates, but other more detailed properties do. In the case of Glauber 
and Metropols dynamics, the low temperature phase can be further divided into two 
regions with different scaling properties of the average trapping time as a function of 
the temperature. Overall, our results rationalize and link the empirical findings about 
correlations between the energy of the minima and their degree, and should stimulate 
further investigations on this issue. 

Our work has also interesting applications in terms of diffusion phenomena on 
complex networks, and shows that non trivial transition rates can lead to a very 
interesting phenomenology. Usual random walks lead to a higher probability for the 
random walker to be in a large degree node (oc kP(k)), with respect to the random choice 
of a node (oc P(k)); here, the models we have studied can lead to various stationary 
probabilities, such as for instance a uniform coverage which does not depend anymore 
on the degree. Interestingly, the biased random walks among traps that we have studied 
can even display a phase transition phenomenon, as either a temperature parameter or 
the network's properties are changed, with the possible presence of a glassy phase with 
slow dynamics. 
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